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, Abstract 

Several binary systems that contain a massive star have been detected in both the radio band and at very 
high energies. In the dense stehar photon field of these sources, gamma-ray absorption and pair creation 
are expected to occur, and the radiation from these pairs may contribute significantly to the observed 
radio emission. We aim at going deeper in the study of the properties, and in particular the morphology, of 
the pair radio emission in gamma-ray binaries. We apply a Monte-Carlo code that computes the creation 
. location, the spatial trajectory and the energy evolution of the pairs produced in the binary system and 

its surroundings. The radio emission produced by these pairs, with its spectral, variability and spatial 
characteristics, is calculated as it would be seen from a certain direction. A generic case is studied first, 
Q^' and then the specific case of LS 5039 is also considered. We find that, confirming previous results, the 

I \ secondary radio emission should appear as an extended radio structure of a few milliarcseconds size. This 

2 ■ radiation would be relatively hard, with fluxes up to ~ 10 mJy. Modulation is expected depending on the 

^ I gamma-ray production luminosity, system eccentricity, and wind ionization fraction, and to a lesser extent 

^ . on the magnetic field structure. In gamma-ray binaries in general, the pairs created due to photon-photon 

interactions can contribute significantly to the core, and generate an extended structure. In the case of 
LS 5039, the secondary radio emission is likely to be a significant fraction of the detected core flux, with 
a marginal extension. 
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1. Introduction et al. 2010; Falcone et al. 2011; Moldon et al. 2011). 

Under the strong ultraviolet photon field of a massive 

High- mass binary systems are a well established class of star, gamma-ray absorption can take place (e.g. Ford 

very high-energy (VHE) gamma-ray emitters. Five VHE 1984; Protheroe & Stanev 1987; Moskalenko & Karacula 

sources^ have been associated so far to binaries that con- 1994; Bednarek 2000; Boetcher & Dermer 2005; Dubus 

tain a massive star: PSR B1259— 63 (Aharonian et al. 2006a; Boetcher 2007; Orellana et al. 2007; Khangulyan 

2005a); LS 5039 (Aharonian et al. 2005a); LS I +61 303 et al. 2008; Reynoso et al. 2008; Sierpowska-Bartosik & 

%^ . (Albert et al. 2006); Cygnus X-1 (Albert et al. 2007); and Torres 2008; Romero et al. 2010), this process leading to 

^/ HESS J0632+057 (Aharonian et al. 2007; Hinton et al. the creation of secondary (electron-positron) pairs with 

2009; Skilton et al. 2009; Falcone et al. 2010; Falcone et energies ^ 10 GeV and peaking around 30 GeV. 

al. 2011; Moldon et al. 2011). Among these five, three Secondary pairs are created under the stellar photon 

show clear modulation of the VHE emission associated field and, as shown below, are likely trapped by the stellar 

to the orbital motion: PSR B1259— 63 (Aharonian et al. wind magnetic field of strength B^. For small enough val- 

2005a), LS 5039 (Aharonian et al. 2006), and LS I +61 303 ues of and high enough values of the secondary particle 

(Albert et al. 2009); thus, the emitter cannot be too energy most of the pair energy is emitted, via inverse 

far from the star in these sources. The short dura- Compton (IC) scattering of stellar photons, as gamma- 

tion of the VHE flare observed in Cygnus X-1 (Albert et rays with energies above the pair-creation threshold. If 

al. 2007) likely implies that this flare is originated not the gamma-ray opacity coefficient r^^ is ^ 1, the electro- 

too far from the star, whereas in HESS J0632+057 the magnetic cascade may significantly affect the distribution 

situation is less clear (see Skilton et al. 2009; Falcone of secondary pairs in the system. However, for this to be 

efficient the initial gamma-ray energy must be well above 

1 ^, • 1 XI u- the pair creation energy threshold, e::^eth~l/e* (i.e. deep 

^ inere is also the new GeV gamma-ray bmary, ^ t'^at \ i ^ 

IFGL J1018.6-5856, which may have also been detected m the Klem-Nishma -KN- regime), where eth, e and are 

at VHE (e.g. Corbet et al. 2011; de Ona Wilhelmi 2010). the threshold, the gamma-ray and the stellar photon en- 
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ergies in rrieC^ units. In this regime, say for energies of 
the secondary pairs ^ 1 TeV and typical stellar tempera- 
tures ~ (3 — 4) X 10^ K, synchrotron cooling suppresses the 
cascade under a magnetic field (Khangulyan et al. 2008) 

5w ^ 10 (L*/10^^ erg s"^) (i?/10^^cm)"^ G, (1) 

where and R are the star luminosity and distance from 
the star center, respectively. Equation 1 has been derived 
from Eq. 1 in Khangulyan et al. (2008), which gives the IC 
cooling timescale in the KN regime with a 10 eV photon 
field as target. We reproduce the formula here for the sake 
of clarity: 

tKN ^ 1.7 X lO^K/lOOerg cm"^)"^ (£;/lTeV)^-^s, (2) 

where is the stellar radiation energy density, and E 
is the electron energy. For > v^Sttu^, the secondary 
bolometric luminosity output will be dominated by syn- 
chrotron radiation. This implies a magnetic field strength: 

B^>3x 10' (L*/10^Vg s-^)^/'(i?/10^'cm)-i G . (3) 

All this shows that even for magnetic fields well below 
equipartition with radiation, synchrotron losses can effec- 
tively prevent development of electromagnetic cascades. 

Bosch- Ramon et al. (2008a) carried out a detailed semi- 
analytical study of the spectral properties of the sec- 
ondary broadband emission for realistic magnetic field val- 
ues (see also Bosch- Ramon et al. 2008b). In that work, it 
was shown that the radio. X-ray, and GeV fluxes of sec- 
ondary emission could be comparable to those observed 
in gamma-ray binaries, and may even dominate the non- 
thermal output below eth in some cases. These authors 
also suggested that the extended radio emission found in 
LS I +61 303 (Dhawan et al. 2006) may come from sec- 
ondary pairs. Bosch-Ramon (2009) also proposed this ori- 
gin for the emission of the marginally resolved radio core 
of LS 5039 (Ribo et al. 2008). 

Although the radio fluxes can be roughly estimated 
from the semi-analytical study of Bosch-Ramon et al. 
(2008a), their approach does not allow a proper descrip- 
tion of the radio spectrum, lightcurve, and morphology, 
because synchrotron radio emitting particles have long 
timescales and are strongly affected by the medium in- 
homogeneity at spatial scales of ~ R. In that work, the 
complex 3-dimensional (3D) structure of the secondary 
trajectories was simplifled attaching the secondary pairs 
to the stellar wind, with the latter being taken as spheri- 
cally symmetric, with constant velocity, and wind particle 
and magnetic fleld energy density ocl/R'^. The magnetic 
fleld was considered fully irregular. Finally, the adiabatic 
cooling of the secondary pairs was not taken into account, 
which is reasonable at binary spatial scales, but not for 
the computation of that radio emission produced outside 
the binary system, affected by wind expansion. Finally, 
the relevance of free-free absorption was not investigated. 

In this work, we present detailed calculations of the evo- 
lution of the radio emitting secondary pairs in gamma-ray 
binaries. The energy and the spatial evolution of the parti- 
cles in the stellar wind have been computed using a Monte- 
Carlo code that accounts for ionization/coulombian col- 
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lision losses, adiabatic, relativistic bremsstrahlung, syn- 
chrotron, and IC cooling, as well as diffusion and ad- 
vection in the wind to follow the trajectories of sec- 
ondary pairs in the source. The code can also deal with 
anisotropic diffusion when an ordered magnetic fleld is 
dominant, although we do not consider this conflguration 
at this stage. A sketch of the the scenario at the scales 
of the binary system is shown in Fig. 1. Once the par- 
ticle energy and spatial distributions have been obtained 
at a particular orbital phase (accounting for the previ- 
ous orbital history of the system), the radio emission in 
a speciflc direction is calculated, which allows the genera- 
tion of directional radio maps. The radio maps have been 
smoothed using a circular Gaussian in order to roughly 
reproduce what an observer would see using a milliarc- 
second (mas) resolution radio interferometer. This is also 
illustrated in Fig. 1. Spectra and lightcurves, and the 
impact of free-free absorption, have been discussed. We 
focus here on a generic case, since our aim is to give a 
general demonstration of the relevance of secondary ra- 
dio emission, but we provide also with an instance of the 
importance of this phenomenon in a real source, applying 
our calculations to LS 5039. Further and more detailed 
studies for speciflc objects of the role of secondary pairs 
for the radio emission will be presented elsewhere. 

2. Monte-Carlo simulations of secondary pairs in 
gamma-ray binaries 

For an appropriate study of the secondary radio emis- 
sion, detailed calculations of the energy and spatial evolu- 
tion of the secondary pairs are required, accounting for the 
star radiation fleld and the properties of the stellar wind. 
Given the difliculties of treating the magnetic fleld struc- 
ture and adiabatic cooling in the wind in a semi-analytical 
approach, we have performed Monte-Carlo simulations in 
which secondaries pairs are injected, move (diffuse and are 
advected), and cool down in a 3D space. The 3D nature of 
the problem comes from the combination of the advection 
in the wind plus the orbital motion (a signiflcant regu- 
lar 5w-field structure would also require a 3D treatment). 
The problem treated by the Monte-Carlo method is sim- 
ilar to that described by the diffusion advection equation 
(e.g. Jones 1990; Blandford & Eichler 1987). The radia- 
tion timescales {t=\E/E\, where E is the electron energy 
loss rate) adapted to the environments of gamma-ray bi- 
naries are: 

Wnc^ 40 (5/100 G)-^E/1 GeV)-^s, (4) 

he Th^50(ii*/300 ergcm-^)-^(£;/l GeV)"^s, (5) 

tbr - 10^(nw/10^°cm-^)-^s, (6) 

for synchrotron, IC (Thomson) and relativistic 
Bremsstrahlung, respectively (see Bosch-Ramon & 
Khangulyan 2009 for the IC cooling rate accounting also 
for the KN regime used in the calculations), with 
being the wind density (see below), the energy density 
of the target photon fleld, and B the ambient magnetic 
fleld. For the adiabatic and ionization cooling timescales, 
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Fig. 1. Left: Sketch of the considered scenario at the scales of the binary system, looking at the system perpen- 
dicularly to the orbital plane. Secondary pairs are created in the vicinity of the gamma-ray emitter, get trapped 
in the stellar wind through slow diffusion, and are advected by it away from the star. Right: Sketch of the ra- 
dio emission produced by secondary pairs as seen by the observer. The radiation emitted by the advected secondary 
pairs, of spiral shape due to wind advection and orbital motion, will appear smoothed by the radio telescope beam. 



we have adopted: 
tad = (3/2) (R/v^oc) (1 - R^/2R)-\which fori? > R^ 
^ 2 X 10"^ {R/3 X 10^^ cm) (vwoo/2 x 10^ cm s"^)"^ s ,(7) 
and 

tion^2x 10^^£;/nw 

^3.4x 10^(£;/1 GeV)(nw/10^cm-^)-is, (8) 

respectively, where v^oo is the wind velocity at infinity, 
and R^ the stellar radius. 

2.1. Injection of secondary pairs in the binary system 

We have calculated the injected spatial and energy dis- 
tribution of secondary pairs for a given primary gamma- 
ray spectrum and a stellar black-body photon field. This is 
done tracking the path of a photon produced in the emit- 
ter location and moving with a certain direction in the 
binary system. The location in which a given gamma ray 
is absorbed is obtained using the (anisotropic) gamma-ray 
absorption differential probability due to photon-photon 
interactions (point-like treatment of the star; see below): 

^ = (1 - cos^) J de*cr^^ (ee* (1 — cosO))n^{e^,R) , (9) 

where 0, and (7 are the length of the path covered 
by the gamma ray, the gamma-ray /stellar photon inter- 
action angle, the stellar photon specific density, and the 
pair-creation cross section (Gould & Schreder 1967), re- 
spectively. This quantity is integrated along the photon 
path Jq dw to obtain the probabity of interaction up to a 
certain location. In this work, the absorption probability 
is computed along the photon trajectory with small steps, 
and the Monte- Carlo algorithm is applied to define the 
interaction point. 

The injection rate of secondary pairs depends strongly 
on the spectrum and angular distribution of the primary 



gamma rays. To mimic the effect of the anisotropic IC 
scattering in the angular distribution of primary gamma 
rays, a gamma-ray direction probability at injection oc (1 — 
cosO^) has been introduced, where 0^ is the angle between 
the gamma-ray path and the direction from the star to 
the gamma-ray emitter, i.e. gamma rays have a higher 
probability to be emitted towards than away from the star. 
This angular dependence would roughly correspond to an 
average one between the Thomson and KN IC regimes 
(see Fig. 6 in Khangulyan et al. 2008). 

In Fig. 2, we show the distribution of injected sec- 
ondary pairs in a gamma-ray binary for three different 
energies. We have focused here on a generic case, in 
which we have adopted a circular binary system (e = 0) 
of one week period, formed by a massive, primary star 
of mass and compact object of mass Mx, such that 
+ Mx = 22 Mq. This yields an orbital separation 
distance of Rorh = 3 x 10^^ cm. The inclination of the 
system has been fixed to i = 45°, and the phase 0.0/0.5 
has been taken at the superior /inferior conjunction of the 
compact object (supc/infc). The main star parameters 
have been fixed to = 6 x 10^^ erg s~^, R^ = 10^^ cm 
and = 3 X 10^ K, where is the star temperature. The 
gamma-ray emitter has been assumed to be point-like, at 
the compact object location, with Rq = i^orb, although the 
actual location and size of the VHE emitter in gamma-ray 
binaries has not been defined yet, being possible to have 
an extended production region not only in the jet but also 
in the colliding wind scenario (see, e.g., Bogovalov et al. 
2008). No eccentricity and a fixed i-value have been taken 
for the sake of simplicity, since our aim is to illustrate the 
main characteristics of the secondary radio emitter. 
2.1.1. Treatment of the star 

In these calculations, we have adopted a point-like ap- 
proximation for the star (e.g. Dubus 2006a; see also 
Sect. 4.2.1 in Bosch-Ramon & Khangulyan 2009). This 
approximation fails close to the stellar surface, at a dis- 
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Fig. 2. Projection in the orbital plane of the spatial distribution of the secondary pairs injected at three different en- 
ergies: E = 10 GeV (left panel), 300 GeV (central panel) and 10 TeV (right panel). This has been computed using 
Monte-Carlo simulations. The gamma-ray emitter is located at {x,y) = (0, 3 x lO"*^^) cm, and the star at (x,y) = (0,0) cm. 
Most of the secondary pairs are produced towards the star, with a plateau in their concentration up to the dis- 
tance at which r^^ ~ 1. This distance changes with energy ^ im^miimi-iim nrr,iinH ff^^^r inn Ctf^Al W mr,ef 
of the interaction angles, corresponding to a peak in the pair 



tance comparable to the radius of the star R^^ i.e. R ^ 
2R^. This region wih be relevant for gamma rays with 
pair-production mean-free paths comparable to the dis- 
tance Re between the star and the gamma-ray emitter. 
Such gamma rays can be absorbed in a volume ~ 47ri?g/3, 
with a sub-region affected by the star finite-size of volume 
^ 7r(2i?*)^ X R^. Thus, even for the most sensitive inter- 
val of gamma-ray energies, the point-like approximation of 
the stellar photon field fails for a fraction of the total rel- 
evant volume ~ 3{R^/Re)^ ~ 0.11 and 0.38 for Rq = 3R^ 
and 2i?*, respectively. Therefore, given a rather broad 
distribution of primary gamma rays and for Rq^2R^^ 
the point-like approximation for the star would have an 
accuracy at least of the order of a 10%. This estimate is 
suported by semi- analytical calculations of the pair injec- 
tion rate density, whose results are presented in Fig. 3. 
The calculation is done computing first the gamma-ray 
flux that reaches a certain point in the binary system. 
Then, the rate density of pairs of 100 GeV (dominant in 
radio) injected in each point is computed, accounting for 
the energy differential cross section of pair creation (see 
Bosch- Ramon et al. 2008a and references therein). This 
calculation have been performed for a generic case and for 
LS 5039 (see Table 2). In Fig. 4, we show the fraction of 
primary gamma-rays absorbed in the system, as a function 
of the gamma-ray energy, for two different distances be- 
tween the gamma-ray emitter and the star: Rq = 3R^ and 
Rq = 2R^^ with the properties of the star chosen as in the 
generic case. As seen from these figures, accounting the 
stellar finite size or adopting a point-like star approxima- 
tion yield almost the same results. This permits to safely 
adopt the point-like approximation for the Monte- Carlo 
calculations presented below, shortening significantly the 
calculation time. 

2.1.2. Primary gamma-ray spectrum and flux 

To perform the calculations corresponding to a generic 
gamma-ray binary, the photon index of the produced 
gamma rays has been taken F = 2.5, and the unab- 
sorbed differential photon rate at 1 TeV is N^eV = 1.4 x 
10^^ TeV~^ s~^, which corresponds to a specific pho- 
ton flux nxev = 3 X 10"^^ TeV"^ s"^ cm"^ at a dis- 




E, TeV 

Fig. 4. Fraction of the primary gamma-rays absorbed in the 
stellar photon field. The gamma-ray emitter was assumed 
to have a (1 — cos ^)-angular dependence (see Sect. 2.1 for 
details). The calculations were performed for the star prop- 
erties given in Table 1. The dotted line corresponds to the 
point-like star approximation and a star-gamma-ray emitter 
distance of Re = 3i?*; the solid line corresponds to the stel- 
lar finite size case and a star-gamma-ray emitter distance of 
Re = 3i?* ; the dash-dot-dot line corresponds to the point-like 
star approximation and a star-gamma-ray emitter distance of 
Re = 2R^] the dashed line corresponds to the stellar finite 
size case and a star-gamma-ray emitter distance of Re = 2R^ . 

tance oi d = 2 kpc. These values are assumed to be or- 
bital phase independent. In Fig. 5, we show the primary 
gamma-ray emission between 20 GeV to 2 TeV in the ob- 
server direction, for the orbital phases 0.0 (supc), 0.25, 
0.5 (infc), and 0.75. The thin lines correspond to the in- 
trinsic spectra, whereas the thick ones take into account 
photon-photon absorption. The computed absorbed spe- 
cific photon fluxes are, at 1 TeV: nxeV ^ 10~^^ (phase 0.0), 
2 X 10-12 (phase 0.25/0.75) and 5 x 10-^^ TeV-^ s-^ cm-^ 
(phase 0.5). 

2.2. Energy and spatial evolution of secondary pairs 

2.2.1. Transport of secondary pairs 

A secondary pair is injected where the gamma ray was 
absorbed, and the subsequent locations are computed ac- 
counting for the local direction and strength of the wind 
magnetic field, which is assumed to have an ordered and 
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Fig. 3. Maps of the injection rate per volume unit for injected pairs of 100 GeV, dominant in radio. This has 
been computed semi-analytically. Two cases are presented: a generic case (top) with the parameter values of 
Table 1, and a case with the properties of LS 5039 around periastron (bottom; see Table 2). The calculation ac- 
counting for the star finite size is shown at the left, and the point-like approximation at the right. Units are 
the logarithm of the rate density normalized to the peak value. The star is located to the left in both maps. 
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Fig. 5. Primary gamma-ray emission between 20 GeV and 
2 TeV in the direction to the observer, for four dif- 
ferent orbital phases: 0.0 (solid line), 0.25/0.75 (long- 
dashed line), and 0.5 (dotted line). The thin lines 
show the unabsorbed spectra, and the thick ones show 
the same but accounting for photon-photon absorption. 



a disordered component. We consider the stellar rotation 
axis to be perpendicular to the orbital plane. 

The ordered wind magnetic field can be described as 
Bw(r,(/)) = ^r^r + ^060, where 

B,^B^{R^/Rf, (10) 

Bcj, ^ B^(v^^/v^oo){R*/R) , (11) 

and B^ and 

are the magnetic field at the star surface 
and the initial azimuthal wind velocity, respectively. For 
the calculations, B^ is fixed to a value of 200 G (see Bosch- 
Ramon et al. 2008b and references therein), high enough 
to suppress the electromagnetic cascade. We neglect here 
the region within the Alfven radius from the star, typically 
confined to a region < 2 i?* , since almost all secondary 
pairs are created farther from the star. The wind velocity 
can be expressed as Vw(r,(/)) = '^wr^r + '^^w^e^, where: 

(i?) (12) 

and the azimuthal one: 
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(13) Table . Binary system properties 
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VwsiR) ~ 0.1 {R,/R). 

The wind density depends on R as 

n^{R) = M^/Att v^r{R) , (14) 

where is the star mass-loss rate. In this work, we have 
taken = 10"^ M© yr"^ and v^oo = 2.5 x 10^ cm s"^ 
The magnetic field and the wind model adopted here are 
based on that presented in Usov & Melrose (1992). 

The irregular component of the magnetic field, SB^ al- 
lows particles to effectively diffuse in the stellar wind. 
The particle mean free paths parallel and perpendicu- 
lar to the magnetic lines are Ay = rjVg and = ^g/^? 
where rg = 3 x 10^(£;/100GeV)(5w/l G)-^ cm is the elec- 
tron gyro radius, and 7^ is > 1 and relates the irregular and 
the wind magnetic field component through B^^ = y/r]6B. 
(thus the ordered component strength is ^/r] — ISB). The 
spatial power spectrum of SB is assumed to be flat, i.e. rj 
does not depend on energy, and the diffusion coefficients 
are proportional to the Bohm one: k.\\^± = A||^^c/3. 

The diffusion approximation is applicable if secondary 
pairs deflect significantly before cooling down, i.e. = 

A||/c <C tcool- 

5w>10-^(^/10^)(^/l GeV)(tcooi/10'cm s-^)-^G.(15) 

In addition, the particle mean free path parallel to the 
magnetic lines (always longer than the perpendicular one) 
should be much smaller than the system typical scale R: 

X\\/R^10-^{E/1 GeV)(5w/lG)-^ 

X {r]/10^)-^{R/3xlO^Y^. (16) 

The simulation time steps have to fulfill as well these con- 
ditions and in addition be ^ for the diffusion approxi- 
mation to be suitable. All this has been accounted for in 
the calculations. 

The diffusion in the stellar wind has been accounted 
in each step j of the secondary evolution as a location 
shift, AX;^.^ = AXyen+AX-^e^, which corresponds to 
the diffusion displacement vector of the particle after a 
time interval dP , accounting for the diffussion coefficients 
Since particles are confined to the stellar wind, they 
will be also advected by the wind. This advection motion 
has been accounted for introducing an additional shift in 
the particle location after each time step j by AX^^^ = 
w^dP . The total spatial shift of the secondary location 
per time step is therefore AX;^-^ + AX^^^. 

For the sake of simplicity, at this stage we have assumed 
a fully disordered magnetic field with 77 = 1. We note that, 
as long as Ay <C R^ i.e. for r] <C 10^ (see Eq. 16), radio 
emitting secondary pairs will basically move anchored to 
the stellar wind. In Table 1, the parameter values are 
given for the generic case studied in this work. 
2.2.2. Secondary particle energy evolution 

Once created at time t*, the secondary pair trajectory 
is calculated by the Monte-Carlo code up to time i.e. 
the time when the energy and spatial distributions of pairs 
and their radiation are computed. The energy evolution 
of secondary pairs is calculated accounting for the en- 
ergy losses along the trajectory taking into account all 



Parameter [units] 


Symbol 


Value 


Stellar radius [cm] 


R^ 


10^2 


Orbit separation distance [cm] 




3 X 10^^ 


Emitter-to-star distance [cm] 


R^ 


3 X 10^^ 


Sitplli^T' tPTYi"nPT"i^ tiiyp T\ 

kJ UCiidl UCiiiJJCl Ct U U.1 c 


T 


S X 1 0^ 


tot?il tl"\7'C;tPTn TYli^Cltl A/fr-^\ 
\j\J\jCX>L oVoUCiii iiidoo iHlQ 


A/f ^ A/f V 


22 


Star surface magnetic field [G] 


B^ 


200 




M 


10~^ 


Wind speed at infinity [cm s~^] 


Voo 


2.5 X 10^ 


Distance [kpc] 


d 


2 


Superior conjunction 


supc 


0.0 


Inferior conjunction 


infc 


0.5 


Inclination angle [°] 


i 


45 


Irregular magnetic field fraction 


V 


1 


1 TeV specific flux [TeV-^ s"^] 


riTeV 


3 X 10-^1 


Gamma-ray photon index 


r 


2.5 



the processes in play: ionization, adiabatic, relativistic 
bremsstrahlung, synchrotron, and IC cooling. Thus, at 
each time step j, a certain amount of energy AE^ = EdP 
is rested to the particle energy at j — 1 , where dP is such 
that AE^ <c£^-^~^, and E includes the cooling mechanisms 
specified above. We note that the energy losses depend 
not only on the particle energy, but also on the particle 
position in the system. 

In the region close to the star {R ~ few i?*), IC cooling 
dominates for basically all the particle energies, although 
at this distance there are none or very few radio emit- 
ting particles. Farther out, adiabatic cooling in the wind 
becomes faster. The energy distribution of radio secon- 
daries almost does not depend on the original gamma-ray 
spectrum, since radio emitting secondary pairs are mostly 
produced around eth and then cool down. Dominant IC 
cooling yields an electron energy distribution below the 
pair creation threshold oc whereas adiabatic cooling 
renders oc E~'^. 

In the Monte-Carlo simulation we have considered an 
equal number of primary gamma-rays (10^ and 2 x 10^ 
for the generic and LS 5039 cases, respectively) per loga- 
rithmic energy bin with the angular distribution described 
in Sect. 2.1. The injection process has been implemented 
during the relevant time interval, starting early enough be- 
fore the time t, to account for all the relevant secondary 
pairs (say, at t^ G [t*min,^))- The orbital motion of the 
primary gamma-ray source has been accounted in the cal- 
culations of the injection rate. The Monte-Carlo code cal- 
culates the secondary pair energy and space distribution 
at the time t keeping record of the time of creation t^ 
and the initial pair energy. To compute the final radiative 
output of the secondary pairs one has to account for the 
adopted gamma-ray injection spectrum (F = 2.5) and the 
corresponding normalization, i.e. one needs to know how 
many particles correspond in reality to one particle of the 
simulation. 
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An analytical code has been used to compute the total 
amount of secondary pairs generated per energy bin for 
the primary gamma-ray distribution dNi^/dj. This value, 
divided by the number of secondary particles per energy 
bin injected in the simulation, gives the weight of an elec- 
tron/positron from the evolved distribution of particles. 
This weight depends on the secondary particle energy, and 
may vary with the orbital phase, i.e. with injection time 
t^ . 

Pairs created at different epochs can contribute to 
the radio emission at the same frequency at the final time 
t, since pairs may have different initial energies and spatial 
trajectories and therefore cooling times differ as well. 

In Figs. 6 and 7, we show the (evolved) secondary loca- 
tion distribution and density map, respectively, projected 
in the observer plane for four different phases, 0.0, 0.25, 
0.5, and 0.75. 

2.3. Secondary pair radiation 

From the energy and the spatial distribution of par- 
ticles at a specific orbital phase, we have computed the 
synchrotron emission convolving the particle energy distri- 
bution with the one-particle power function. Computing 
the fluxes from 2D small regions or bins in the observer 
plane (image pixels), i.e. integrating in depth, permits 
to extract radio morphologies. The radio emission has 
been calculated in the optically thin approximation of syn- 
chrotron radiation, valid at frequencies above few GHz for 
the relevant spatial scales. For a 10 mJy source with the 
self- absorption turnover at ~ 1 GHz, the radio emitter 
size would be constrained to ^ 5 x 10^^(5/1 G)^/^ cm 
(Bosch-Ramon 2009) (where B would be the radio emit- 
ter characteristic magnetic field), or ~ 1 mas at few kpc. 

Typical radio fluxes of about 20 mJy are obtained at 

5 GHz for the parameter values adopted here. Note that, 
fixing F, the secondary radio flux is oc N^^y /d^ . The spec- 
tral shape of the specific flux is almost flat, Fj^ oc constant, 
since adiabatic losses dominate the secondary pair evolu- 
tion. Part of the radio emission originates nevertheless in 
regions in which IC cooling dominates, which softens the 
spectrum. Free- free absorption, of coefficient rff^^, could 
be relevant. Adopting a simple prescription for rff^^ which 
does not depend on the orbital phase in a circular orbit 
(see Eq. 1 in Bosch- Ramon 2009, adapted from Rybicki 

6 Lightman 1979), we find that for ionization fractions 
-^ion ^0.1, the radio fluxes should be already significantly 
affected. Another effect of the free-free absorption is a 
strong spectral hardening, although for Xion ^ 1 a com- 
bination of the energy and the spatial distribution of sec- 
ondaries yields again a flat spectrum. All this is shown 
in Fig. 8. Regarding the magnetic field geometry, fixing 
77 = 1 implies that the synchrotron power does not depend 
on the magnetic field orientation and therefore the orbital 
phase. For 7^ > 1, changes in the B^-Yme of sight angle 
could introduce moderate orbital flux variability. For a 
dominantly toroidal magnetic field, synchrotron flux vari- 
ations would be of the order of sin (7r/2)/sin^(i) (2 for 
z = 45°), with two dips and two enhancements along the 
orbit when the By^-\m.e of sight angle reached i and 7r/2, 



- x.=o 
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Fig. 8. Computed spectra of the whole radio emit- 
ter for three ionization fractions, -^ion = (sohd 
hne), 0.1 (long-dashed line) and 1 (dotted line). 

respectively. 

The radiation image can be computed as it would be 
seen by the observer projecting the 3D emitter in the ob- 
server plane and dividing then in 2D bins. From that, syn- 
thetic radio maps can be produced convolving the fluxes 
in these bins with a circular Gaussian of full width half 
maximum (FWHM) of 1 mas. The result will be similar to 
an image obtained by a radio interferometer with angular 
resolution of 1 mas. In Figs. 9, 10 and 11, we show three 
sets of 5 GHz radio images at four different phases: 0.0 
(supc), 0.25, 0.5 (infc), and 0.75. Fig. 9 corresponds to the 
raw image, i.e. without Gaussian convolution. There is 
an expectable similarity in shape between the radio sec- 
ondary (Fig. 7) and the flux spatial distribution. The 
image convolved with a Gaussian with FWHM=1 mas 
is presented in Fig. 10, and Fig. 11 shows the residuals 
when extracting a point-like source (the FWHM=1 mas 
Gaussian) to the convolved image both with the same flux. 
Typical noise levels for such an interferometer could be 

0.05 — 0.1 mJy/beam, slightly below the predicted ex- 
cesses (see the intensity scale in Fig. 11). The impact of 
free- free absorption has not been included in these plots, 
which would consist on a reduction of the fluxes in the 
image core but keeping an excess in the periphery. 

We have focused in this work on the synchrotron radio 
emission. We note nevertheless that for the parameter 
values adopted here, the synchrotron X-ray and IC GeV 
fluxes would be both about few 10~^^ erg s~^ cm~^. A 
more detailed study of the higher energy secondary radi- 
ation will be presented elsewhere. 

3. Application to LS 5039 

We have applied our Monte- Carlo code and radiation 
model presented above to LS 5039. LS 5039 is a high-mass 
X-ray binary, detected in VHE gamma rays (Aharonian et 
al. 2005b; Aharonian et al. 2006), which has been observed 
in radio with different VLBI instruments for more than a 
decade (Marti et al. 1998; Paredes et al. 2000; Paredes 
et al. 2002; Ribo et al. 2008). The work by Ribo et al. 
2008 presents VLBA observations at 5 GHz in which the 
radio core appears as marginally resolved, with two small 
dipolar elongated extensions that show a moderate change 
in the position angle. The core and extension fluxes are 
18-20 mJy and 6-4 mJy, respectively, having the whole 
structure an angular size of several mas. The radio struc- 
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Fig. 6. Computed spatial distribution of the secondary pairs, projected in the observer plane (i = 45°), for four differ- 
ent phases: 0.0 (top, left), 0.25 (top, right), 0.5 (bottom, left), and 0.75 (bottom, right). About 10^ particles have 
been injected (black spots), among which about 10"^ particles have GHz synchrotron emitting energies (light blue spots). 
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Fig. 7. Density map of the radio secondary populations shown in Fig. 6. Intensity units are arbitrary. 
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Fig. 9. Computed image of the 5 GHz radio emission, in the direction to the observer, in units of mJy per beam. 
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Fig. 10. Gaussian convolved image of the 5 GHz radio emission shown in Fig. 9. 
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Fig. 11. Residuals of the images presented in Fig. 10 after subtracting a point-like (Gaussian) source with the same flux. 



tures found in the source have been associated either to a 
jet (Paredes et al. 2000) or a shocked pulsar wind (Dubus 
2006b). The binary system consists of an 06.5 main se- 
quence star and a compact object of unclear nature (neu- 
tron star or black hole; see Casares et al. 2005). The 
system semi-major axis is a = 2.2 x 10^^ cm, with eccen- 
tricity e = 0.34 and orbital period P = 3.9 days. The phase 
0.0 corresponds to the periastron passage, and phases 0.05 
and 0.67 to supc and infc of the compact object, respec- 
tively. The inclination angle of the orbit is not known 
and could be in a range between i ^ 15° — 75°. The stellar 
mass- loss rate has been taken to be ^ 3x 10~^ Mq/jt. 
The ephemeris, orbital parameters, system-observer ge- 
ometry, star and stellar wind properties, and inclination 
angle constraints have been taken from Casares et al. 2005 
and Aragona et al. 2009 (see also Sarty et al. 2010). The 
system orientation in the plane of the sky, unconstrained 
from observations, has been taken as in Fig. 6 of Aragona 
et al. 2009. Following the approach for the generic case, 
i has been fixed to 45° (though our results should not 
change qualitatively for i values 45dzl5°). The values 
of and rj have been taken also as in the generic case, 
200 G and 1, respectively. The parameters adopted for 
LS 5039 are presented in Table 2. 

Given the uncertainties in the source properties, and the 
known complexities of the orbit and VHE fiux/photon- 
index light curve, we do not aim in this section at carrying 
a thorough analysis of the different possibilities, what is 
left for future work. Otherwise, we intend to give two 



possible instances of the source behavior, and for that 
we focus on two phases, 0.2, slightly after supc, and 0.8, 
slightly after infc. The distribution of the gamma rays in- 
jected in the system has been chosen such that the source 
should show absorbed fiuxes and photon indices similar to 
the observed ones (Aharonian et al. 2006) in the direction 
to the observer. For this, we have deabsorbed the ob- 
served VHE spectra (see, e.g., Boetcher 2007) assuming a 
point-like emitter in a certain location, as explained below, 
emitting in the energy range from 20 GeV up to 2 TeV. 
Since the observed lightcurve shows a relatively smooth 
behavior around phases 0.45-0.9 and 0.9-0.45 (Aharonian 
et al. 2006), for simplicity we have adopted typical specific 
fiuxes and photon indices to compute the radio emission: 
riTeV = 10"^^ TeV"^ s"^ cm-^ and F = 2.5 (phase 0.2), 
and 3 x 10"^^ TeV"^ s"^ cm-^ and F = 2 (phase 0.8). 
Since the location of the emitter is not well constrained 
(Khangulyan et al. 2008; Bosch- Ramon et al. 2008b), two 
different heights have been chosen that optimize the radio 
fiuxes without entering in strong confiict with the X-ray 
fiuxes. At phase 0.2, we have located the emitter at a 
height of 10^^ cm, and at 2 x 10^^ cm at phase 0.8, perpen- 
dicularly to the orbital plane, in the observer side of that 
plane. There is no additional emitter assumed at a loca- 
tion symmetric with respect to the orbital plane, although 
if gamma-rays were produced in a symmetric bipolar jet, 
this component should be also included in the calcula- 
tions. However, given the expected higher wind free- free 
opacity for this component, its contribution is expected 
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Table . LS 5039 properties 



Parameter [units] Symbol Value 

Stellar radius [cm] 7 x 10^^ 

Orbital semi-major axis [cm] Rorh 2.2 x 10^^ 

eccentricity e 0.34 

Stellar temperature [K] 3.8 x 10"^ 
total system mass [M©] + Mx 25 

Star surface magnetic field [G] 200 

Mass-loss rate [M© yr"^] M 3 x 10"'^ 

Wind speed at infinity [cm s~^] Vcx) 2.4 x 10^ 

Distance [kpc] d 2.5 

Superior conjunction supc 0.05 

Inferior conjunction infc 0.67 

Inclination angle [°] i 45 

Irregular magnetic field fraction rj 1 



could have a non-negligible component coming from the 
secondary pairs created in the vicinity of the gamma-ray 
emitter. Actually, the following simple estimate of the 
radio flux already indicates this. Assuming that secondary 
pairs cool mainly through adiabatic cooling in the region 
in which they produce radio emission, which is the case 
for a wide range of parameter values, and a wind velocity 
of 2 X 10^ cm s~^, the expected fluxes can be derived from 
the formula: 

F3 GHz ^ {l/47r(P)L^± (£^radio/^0) (WAsync)rioO (17) 

= 20(LvHE/10''erg s'^) {B^^/lOG) 

x(L*/10^^erg s-^){Re/3 x lO^^cm)"^ ((i/2kpc)- Vjy , 

where £^radio is the energy of the radio emitting parti- 
cles, Tioo the angle averaged opacity at gamma-ray en- 
ergies around 100 GeV, tgync the synchrotron cooling 
timescale, Eq ~ 30 GeV the secondary typical injection 
energy, Lvhe the VHE luminosity, = B(f){RQ)^ and 
d the distance to the source. In Eq. (18), the prod- 
uct (L*/10^^erg s-^) (i?e/3 x lO^^cm)"^ should be substi- 
tuted by ~ 2 for rioo > 1 (and rioo fixed to 1). In sources 
with parameters such that F5 ghz ^ 1 mJy, secondary ra- 
dio emission should not be neglected. The most crucial 
point here is whether magnetic fields of the order of 10 G 
can be found at few R^ from the star. 

As shown in this work, secondary radio emission from 
gamma-ray binaries would not be only relevant in flux, 
but they could also present resolvable extensions. In gen- 
eral, the spectra can be flat or softer depending on 
(adiabatic vs IC cooling dominance), and the absorbed 
gamma-ray luminosity and the radio emitter conditions, 
in particular and B^^ along the orbit. In the par- 
ticular case of LS 5039, we find that a dominant frac- 
tion of the marginally extended core of LS 5039 (Ribo et 
al. 2008) could be of secondary origin, with a changing 
relatively soft spectrum. Since the computed additional 
extension appears marginal, the observed elongated emis- 
sion at 5 GHz beyond the core in LS 5039 would be likely 
related to an intrinsic primary radio emitter. 

Regarding morphology, we find that a spiral-like shape 
would be expected, although typically the small angular 
size of the radio source presently prevents detailed imag- 
ing. Such a spiral-like shape was also predicted in Bosch- 
Ramon et al. (2008a), although the treatment there was 
simpler than in the present work. It is worthy noting that 
the wind conditions could differ from those adopted here, 
and then the shape, and even the radio emitter extension, 
may change. More collimated radio structures, or a mag- 
netic field increase downstream the wind leading to more 
extended radio emission, cannot be discarded. 

It is remarkable that the spatial distribution of the radio 
emitting secondary pairs is strongly determined by their 
injection, which depends on the stellar photon field and 
the gamma-ray emitter location, and by the stellar wind 
inhomogeneities and anisotropics, but not by B^. Only a 
very regular B^^ with very large /^-values, could modify 
significantly the distribution of radio emitting secondary 
pairs, since then these particles would move at c along 



to be smaller than the one considered here. We have 
computed the secondary injection and evolution starting 
a phase interval of 0.5 earlier than the phases in which 
the radiation is evaluated, i.e. 0.7 0.2 and 0.3 0.8. 
However, most of the radio emission, including the out- 
ermost region of the emitter, comes from secondary pairs 
injected within a phase interval of 0.2 earlier, due to adi- 
abatic cooling. 

The resulting images of the raw radio emission, of the 
source convolved with a Gaussian, and the residuals af- 
ter extracting a point-like source are presented in Fig. 12. 
The figure shows that the radiation, with a typical emit- 
ter size of ^ 1 mas, would be difficult to resolve, although 
the bigger size of the system around infc may lead to a 
marginal extension. The obtained total fluxes at 5 GHz 
are about 6 and 10 mJy at phases 0.2 and 0.8, respec- 
tively. The complex flux evolution along the orbit of 
the primary (deabsorbed) gamma-ray spectra yield dif- 
ferent spectra from those of the generic case, obtaining 
Fj, oc (0.2) and oc (0.8). This is due to secondary 
pairs injected at different phases, with different absorbed 
luminosities, produce the radiation at different frequencies 
at the same orbital phase. The computed luminosities in 
the range 1-10 keV are of about 10^^ erg s~^, similar to 
those observed in this energy range (e.g. Bosch-Ramon et 
al. 2007; Takahashi et al. 2009). These X-ray fluxes al- 
ready constrain the fluxes of the radio emission to several 
mJy at most. These results are nevertheless sensitive to 
the radial profile of B^^ which is not well known, thus we 
do not treat here the high-energy output of the secondary 
emission. The impact of free-free absorption in the wind 
has not been considered either, but this process could re- 
duce significantly the secondary radio emission mainly for 
phases around supc under high ionization fractions. 

4. Discussion 

This work shows that the radio emission detected from 
compact gamma-ray binaries that contain a massive star 
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Fig. 12. Computed raw (top), Gaussian convolved (middle), and residual images (bottom) of the emis- 
sion at 5 GHz from LS 5039 in the direction to the observer for the phases 0.2 (right) and 0.8 (left). 
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the dominantly toroidal B^, i.e. leaving the system with 
a speed c5wr/^w(/), between and c. For reasonable 
and r^-values, particles are trapped in the wind, and 
their energy distribution is determined by IC and adia- 
batic cooling. The magnetic field strength is then rele- 
vant only for the synchrotron fluxes (recall F oc ^w), and 
its level of (ir) regularity, 77, for variability. Polarization 
would also be affected by 7^, and polarization observations 
may probe the B^ geometry. However, the high density of 
cold free electrons in the stellar wind could induce strong 
Faraday rotation. A proper study of the impact of this on 
the final polarization degree and angle requires a devoted 
investigation. 

System eccentricity, a significant regular magnetic field 
component, and free- free absorption would be sources of 
radio variability and should be considered when study- 
ing specific cases. Changes in the gamma-ray luminosity 
would also lead to a smooth modulation of the secondary 
radio emission. To finish with, it seems likely that in 
gamma-ray binaries, the radio emission from inner regions 
can be, if not dominated, significantly contaminated by 
secondary radiation. Secondary pairs should not be ne- 
glected when understanding the broadband non-thermal 
emission in compact gamma-ray sources, nor the physical 
connection of the radiation produced in different bands. 
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